import os, math
import numpy as np
import pandas as pd
import plotly.express as px
import plotly.offline as pyo
import matplotlib.pyplot as plt
%matplotlib inline
import plotly.io as pio
pio.renderers.default = "notebook+pdf"

data_dict = {}
header = {}
i = 0
for filename in os.listdir():
    if "_data" in filename:
        if "_spect" in filename:
            data_dict[filename[:-15]] = pd.read_csv(filename)
        else:
            data_dict[filename[:-9]] = pd.read_csv(filename)
            
    elif "_header" in filename:
        head = pd.read_csv(filename, names = ["Name", "Value"])
        if "_spect" in filename:
            header[filename[:-17]] = head.at[16,"Value"]
        else:
            header[filename[:-11]] = head.at[16,"Value"]
        i += 1
        
print("Example of dataframe:")
data_dict["background"]

Example of dataframe:

	Channel	Energy	Counts
0	0	-62.744125	0
1	1	-60.751958	0
2	2	-58.759791	0
3	3	-56.767624	0
4	4	-54.775457	0
...	...	...	...
1019	1019	1967.274151	3
1020	1020	1969.266319	1
1021	1021	1971.258486	4
1022	1022	1973.250653	1
1023	1023	1975.242820	1

1024 rows × 3 columns

import plotly.graph_objects as go
fig3 = go.Figure()
fig3.add_trace(go.Scatter(x=data_dict["180s_2017_sample6"]["Energy"],y=data_dict["180s_2017_sample6"]["Counts"],
                    mode='lines',
                    name="180s_2017_sample6"))
fig3.add_trace(go.Scatter(x=data_dict["BlockT_180s_2017_sample6"]["Energy"],y=data_dict["BlockT_180s_2017_sample6"]["Counts"],
                    mode='lines',
                    name='BlockT_180s_2017_sample6'))
fig3.add_trace(go.Scatter(x=data_dict["background"]["Energy"],y=data_dict["background"]["Counts"]*float(header["180s_2017_sample6"])/float(header["background"]),
                    mode='lines',
                    name='Background'))
#fig3.update_yaxes(type="log")

def ECalib(ener):
    return (ener+62.74412532637075)/1.9921671018276763

def integPeak(xmin,xmax,df):
    minIn = round(ECalib(xmin))
    maxIn = round(ECalib(xmax))
    return df["Counts"][minIn:maxIn].sum()

Radioactive Decay

def halflife(N0, Nt, t):
    return t/math.log(Nt/N0,0.5)

Samples 5_2013 and 7_2017, Co60, theoretical HL of 5.271y

fig3 = go.Figure()
fig3.add_trace(go.Scatter(x=data_dict["2013_sample5"]["Energy"],y=data_dict["2013_sample5"]["Counts"]*float(header["2017_sample7"])/float(header["2013_sample5"]),
                    mode='lines',
                    name='2013_sample5'))
fig3.add_trace(go.Scatter(x=data_dict["2017_sample7"]["Energy"],y=data_dict["2017_sample7"]["Counts"],
                    mode='lines',
                    name='2017_sample7'))
fig3.update_layout(xaxis_title='Energy [KeV]',
                   yaxis_title='Counts')

p1 = integPeak(1120.603,1246.11,data_dict["2013_sample5"]*float(header["2017_sample7"])/float(header["2013_sample5"]))
p1p = integPeak(1120.603,1246.11,data_dict["2017_sample7"])
print("Integral of peak 1: "+str(p1))
print("Integral of peak 1': "+str(p1p))

p2 = integPeak(1256.07,1437.58,data_dict["2013_sample5"]*float(header["2017_sample7"])/float(header["2013_sample5"]))
p2p = integPeak(1256.07,1437.58,data_dict["2017_sample7"])
print("Integral of peak 2: "+str(p2))
print("Integral of peak 2': "+str(p2p))
print()

print("HL from peak (-1): " + str(halflife(p1p, p1, 3)))
print("HL from peak (0): " + str(halflife(p1p, p1, 4)))
print("HL from peak (+1): " + str(halflife(p1p, p1, 5)))
print()

print("HL from peak (-1): " + str(halflife(p2p, p2, 3)))
print("HL from peak (0): " + str(halflife(p2p, p2, 4)))
print("HL from peak (+1): " + str(halflife(p2p, p2, 5)))

Integral of peak 1: 16623.9933101173
Integral of peak 1': 30989
Integral of peak 2: 13022.771535923754
Integral of peak 2': 23920

HL from peak (-1): 3.338938216649393
HL from peak (0): 4.451917622199191
HL from peak (+1): 5.564897027748988

HL from peak (-1): 3.4200472219342024
HL from peak (0): 4.560062962578936
HL from peak (+1): 5.700078703223671

Samples 6_2013 and 8_2017, Mn54, theoretical HL of 312.2d

s2013 = "2013_sample6"
s2017 = "2017_sample8"

fig3 = go.Figure()
fig3.add_trace(go.Scatter(x=data_dict[s2013]["Energy"],y=data_dict[s2013]["Counts"]*float(header[s2017])/float(header[s2013]),
                    mode='lines',
                    name=s2013))
fig3.add_trace(go.Scatter(x=data_dict[s2017]["Energy"],y=data_dict[s2017]["Counts"],
                    mode='lines',
                    name=s2017))

fig3.update_layout(xaxis_title='Energy [KeV]',
                   yaxis_title='Counts')
#title='Measure of {} and {} Gamma Ray Production, Counts vs. Energy'.format(s2013,s2017),

xmin = 765.9974
xmax = 925.3708
p1 = integPeak(xmin,xmax,data_dict["2013_sample6"]*float(header["2017_sample8"])/float(header["2013_sample6"]))
p1p = integPeak(xmin,xmax,data_dict["2017_sample8"])
print("Integral of peak 1: "+str(p1))
print("Integral of peak 1': "+str(p1p)+"\n")

print("HL from peak (-1): " + str(halflife(p1p, p1, 3)))
print("HL from peak (0): " + str(halflife(p1p, p1, 4)))
print("HL from peak (+1): " + str(halflife(p1p, p1, 5)))

Integral of peak 1: 1711.284378109453
Integral of peak 1': 9781

HL from peak (-1): 1.1928892120792631
HL from peak (0): 1.5905189494390177
HL from peak (+1): 1.988148686798772

print("HL from p1: " + str(halflife(276, 35.12, 4)))

HL from p1: 1.3448526611739258

Samples 7_2013 and 1_2017, Ba133, theoretical HL of 10.53y

s2013 = "2013_sample7"
s2017 = "2017_sample1"

fig3 = go.Figure()
fig3.add_trace(go.Scatter(x=data_dict[s2013]["Energy"],y=data_dict[s2013]["Counts"]*(float(header[s2017])/float(header[s2013])),
                    mode='lines',
                    name=s2013))
fig3.add_trace(go.Scatter(x=data_dict[s2017]["Energy"],y=data_dict[s2017]["Counts"],
                    mode='lines',
                    name=s2017))
fig3.update_yaxes(type="log")
fig3.update_layout(xaxis_title='Energy [KeV]',
                   yaxis_title='Counts')
fig3.show()

######################## Taking the integrals of the peak and calculating Half-Life #########################

xmin = 222.1358
xmax = 423.3446
p1 = integPeak(xmin,xmax,data_dict[s2013]*(float(header[s2017])/float(header[s2013])))
p1p = integPeak(xmin,xmax,data_dict[s2017])
print("First integral: "+str(p1))
print("Second integral: "+str(p1p)+"\n")
print("HL from peak (-1): " + str(halflife(p1p, p1, 3)))
print("HL from peak (0): " + str(halflife(p1p, p1, 4)))
print("HL from peak (+1): " + str(halflife(p1p, p1, 5)))

First integral: 197111.41901922962
Second integral: 299631

HL from peak (-1): 4.965444340086574
HL from peak (0): 6.620592453448765
HL from peak (+1): 8.275740566810956

#Half life from single point
print("HL from p1: " + str(halflife(7604, 4958.303, 4)))

HL from p1: 6.4839059410690885

Samples 8_2013 and 5_2017, Zn65, theoretical HL of 243d

s2013 = "2013_sample8"
s2017 = "2017_sample5"

fig3 = go.Figure()
fig3.add_trace(go.Scatter(x=data_dict[s2013]["Energy"],y=data_dict[s2013]["Counts"]*(float(header[s2017])/float(header[s2013])),
                    mode='lines',
                    name=s2013))
fig3.add_trace(go.Scatter(x=data_dict[s2017]["Energy"],y=data_dict[s2017]["Counts"],
                    mode='lines',
                    name=s2017))
#fig3.update_yaxes(type="log")
fig3.update_layout(xaxis_title='Energy [KeV]',
                   yaxis_title='Counts')
fig3.show()

######################## Taking the integrals of the peak and calculating Half-Life #########################

xmin = 1054.862
xmax = 1162.439
p1 = integPeak(xmin,xmax,data_dict[s2013]*float(header[s2017])/float(header[s2013]))
p1p = integPeak(xmin,xmax,data_dict[s2017])
print("First integral: "+str(p1))
print("Second integral: "+str(p1p)+"\n")
print("HL from peak (-1): " + str(halflife(p1p, p1, 3)))
print("HL from peak (0): " + str(halflife(p1p, p1, 4)))
print("HL from peak (+1): " + str(halflife(p1p, p1, 5)))

First integral: 516.3999696140991
Second integral: 2055

HL from peak (-1): 1.5055875568362662
HL from peak (0): 2.007450075781688
HL from peak (+1): 2.50931259472711

#Half life from single point
print("HL from p1: " + str(halflife(62, 6.54, 4)))

HL from p1: 1.232701471488777

Samples 4_2013 and 4_2017, Co57, theoretical HL of 271.8d

s2013 = "2013_sample4"
s2017 = "2017_sample4"

fig3 = go.Figure()
fig3.add_trace(go.Scatter(x=data_dict[s2013]["Energy"],y=data_dict[s2013]["Counts"]*(float(header[s2017])/float(header[s2013])),
                    mode='lines',
                    name=s2013))
fig3.add_trace(go.Scatter(x=data_dict[s2017]["Energy"],y=data_dict[s2017]["Counts"],
                    mode='lines',
                    name=s2017))
fig3.update_yaxes(type="log")
fig3.update_layout(xaxis_title='Energy [KeV]',
                   yaxis_title='Counts')
fig3.show()

######################## Taking the integrals of the peak and calculating Half-Life #########################

xmin = 56.7859
xmax = 118.5431
p1 = integPeak(xmin,xmax,data_dict[s2013]*float(header[s2017])/float(header[s2013]))
p1p = integPeak(xmin,xmax,data_dict[s2017])
print("First integral: "+str(p1))
print("Second integral: "+str(p1p)+"\n")
print("HL from peak (-1): " + str(halflife(p1p, p1, 3)))
print("HL from peak (0): " + str(halflife(p1p, p1, 4)))
print("HL from peak (+1): " + str(halflife(p1p, p1, 5)))

First integral: 6595.638372747748
Second integral: 54951

HL from peak (-1): 0.9808532444341904
HL from peak (0): 1.3078043259122538
HL from peak (+1): 1.6347554073903172

Radiation Shielding

Aluminium Shield

shields_Al = ["G","J","M","O"]
integrAl = {}
thicknAl = {"G":0.508,"J":1.016,"M":2.032,"O":2.54}
surfDensAl = {"G":141,"J":258,"M":522,"O":655}

fig = go.Figure()

for shield in shields_Al:
    fig.add_trace(go.Scatter(x=data_dict["Block"+shield+"_130s_2017_sample3"]["Energy"],y=data_dict["Block"+shield+"_130s_2017_sample3"]["Counts"],
                    mode='lines',
                    name="Block"+shield))
    integrAl[shield] = integPeak(580.7258,760.0209,data_dict["Block"+shield+"_130s_2017_sample3"])
    
fig.update_layout(xaxis_title='Energy [KeV]',
                   yaxis_title='Counts')
fig.show()

print(integrAl)

{'G': 9984, 'J': 10059, 'M': 9685, 'O': 9765}

Lead Shield

shields_Pb = ["Q","R","S","T"]
integrPb = {}
thicknPb = {"Q":0.8128,"R":1.6256,"S":3.81,"T":6.35}
surfDensPb = {"Q":1120,"R":2066,"S":3448,"T":7367}

fig = go.Figure()
for shield in shields_Pb:
    fig.add_trace(go.Scatter(x=data_dict["Block"+shield+"_130s_2017_sample3"]["Energy"],y=data_dict["Block"+shield+"_130s_2017_sample3"]["Counts"],
                    mode='lines',
                    name="Block"+shield))
    integrPb[shield] = integPeak(580.7258,760.0209,data_dict["Block"+shield+"_130s_2017_sample3"])

fig.update_layout(xaxis_title='Energy [KeV]',
                   yaxis_title='Counts')
fig.show()

print(integrPb)

{'Q': 9144, 'R': 8521, 'S': 7131, 'T': 4923}

integr = {**integrAl,**integrPb}
thickn = {**thicknAl,**thicknPb}
surfDens = {**surfDensAl,**surfDensPb}

testArr = {'Surface Density [g/cm]': [], 'Thickness [mm]': [],'Counts in Peak': []}
a = []
for letter in shields_Al+shields_Pb:
    testArr['Surface Density [g/cm]'].append(surfDens[letter])
    testArr['Thickness [mm]'].append(thickn[letter])
    testArr['Counts in Peak'].append(integr[letter])
    
shieldDf = pd.DataFrame(testArr)
shieldDf

	Surface Density [g/cm]	Thickness [mm]	Counts in Peak
0	141	0.5080	9984
1	258	1.0160	10059
2	522	2.0320	9685
3	655	2.5400	9765
4	1120	0.8128	9144
5	2066	1.6256	8521
6	3448	3.8100	7131
7	7367	6.3500	4923

from plotly.subplots import make_subplots
import seaborn as sns

fig3 = make_subplots(rows=1, cols=2, subplot_titles=("Default Scale", "Logarithmic Scale"))

fig3.add_trace(go.Scatter(x=shieldDf['Surface Density [g/cm]'],y=shieldDf['Counts in Peak'],
                    mode='markers',
                    name="180s_2017_sample6"),col=1,row=1)
fig3.add_trace(go.Scatter(x=shieldDf['Thickness [mm]'],y=shieldDf['Counts in Peak'],
                    mode='markers',
                    name='BlockT_180s_2017_sample6'),col=2,row=1)

fig = px.scatter(shieldDf, x='Surface Density [g/cm]', y='Counts in Peak', trendline="ols")
fig.show()

results = px.get_trendline_results(fig)
results.px_fit_results.iloc[0].summary()

/home/jonas/anaconda3/lib/python3.7/site-packages/scipy/stats/stats.py:1535: UserWarning:

kurtosistest only valid for n>=20 ... continuing anyway, n=8

OLS Regression Results

Dep. Variable:	y	R-squared:	0.986
Model:	OLS	Adj. R-squared:	0.984
Method:	Least Squares	F-statistic:	436.1
Date:	Wed, 18 Nov 2020	Prob (F-statistic):	7.85e-07
Time:	00:08:52	Log-Likelihood:	-53.558
No. Observations:	8	AIC:	111.1
Df Residuals:	6	BIC:	111.3
Df Model:	1
Covariance Type:	nonrobust

	coef	std err	t	P>|t|	[0.025	0.975]
const	1.007e+04	104.686	96.152	0.000	9809.595	1.03e+04
x1	-0.7263	0.035	-20.883	0.000	-0.811	-0.641

Omnibus:	4.563	Durbin-Watson:	2.390
Prob(Omnibus):	0.102	Jarque-Bera (JB):	1.403
Skew:	-1.016	Prob(JB):	0.496
Kurtosis:	3.277	Cond. No.	3.95e+03

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 3.95e+03. This might indicate that there are
strong multicollinearity or other numerical problems.

from plotly.subplots import make_subplots
fig3 = go.Figure()

fig3.add_trace(go.Scatter(x=data_dict["180s_2017_sample6"]["Energy"],y=data_dict["180s_2017_sample6"]["Counts"]-data_dict["background"]["Counts"]*float(header["180s_2017_sample6"])/float(header["background"]),
                    mode='lines',
                    name="180s_2017_sample6"))
fig3.add_trace(go.Scatter(x=data_dict["BlockT_180s_2017_sample6"]["Energy"],y=data_dict["BlockT_180s_2017_sample6"]["Counts"]-data_dict["background"]["Counts"]*float(header["180s_2017_sample6"])/float(header["background"]),
                    mode='lines',
                    name='BlockT_180s_2017_sample6'))
fig3.add_trace(go.Scatter(x=data_dict["background"]["Energy"],y=data_dict["background"]["Counts"]*float(header["180s_2017_sample6"])/float(header["background"]),
                    mode='lines',
                    name='Background'))
fig3.update_yaxes(type="log")